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INTRODUCTION 

The objective of the present research is to develop a general mathemati- 
cal model and solution methodologies for analyzing structural response of 
thin, metallic shell-type structures under large transient, cyclic or static 
thermomechanical loads. Among the system responses, which are associated 
with these load conditions, are thermal buckling, creep buckling and ratchet- 
ting. Thus, geometric as well as material -type nonlinearities (of high 
order) can be anticipated and must be considered in the development of the 
mathematical model. Furthermore, this must also be accommodated in the 
solution procedures. 




SUMMARY OF PROGRESS 

The progress update has been elaborated upon in an interim scientific 
report submitted to the sponsor during the summer of 1986. A second report, 
describing the developed finite element, the computer code, and several 
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applications to cylindrical and spherical shell configurations, is being 
prepared. It is expected that this second scientific report will be sub- 
mitted to the sponsor during the summer of 1987. 

A complete true ab-inito rate theory of kinematics and kinetics for 
continuum and curved thin structures, without any restriction on the magni- 
tude of the strains or the deformation, was formulated. The time dependence 
and large strain behavior are incorporated through the introduction of the 
time rates of the metric and curvature in two coordinate systems; a fixed 
(spatial) one and a convected (material) coordinate system. The relations 
between the time derivative and the covariant derivative (gradient) have been 
developed for curved space and motion, so that the velocity components supply 
the connection between the equations of motion and the time rate of change of 
the metric and curvature tensors. 

The metric tensor (time rate of change) in the convected material 
coordinate system is linearly decomposed into elastic and plastic parts. In 
this formulation, a yield function is assumed, which is dependent on the rate 
of change of stress, metric, temperature, and a set of internal variables. 
Moreover, a hypoelastic law was chosen to describe the thermoelastic part of 
the deformation. 

A time and temperature dependent viscoplastic model was formulated in 
this convected material system to account for finite strains and rotations. 
The history and temperature dependence were incorporated through the intro- 
duction of internal variables. The choice of these variables, as well as 
their evolution, was motivated by phenomenological thermodynamic considera- 
tions. 

The nonisothermal elastic-viscoplastic deformation process was described 
completely by "thermodynamic state" equations. Most investigators (in the 
area of viscoplasticity) employ plastic strains as state variables. Our 
study shows that, in general, use of plastic strains as state variables may 
lead to inconsistencies with regard to thermodynamic considerations. 
Furthermore, the approach and formulation employed by all previous investiga- 
tors lead to the condition that all plastic work is completely dissipated. 
This, however, is in contradiction with experimental evidence, from which it 
emerges that part of the plastic work is used for producing residual stresses 
in the lattice, which, when phenomenologically considered, causes hardening. 



Both limitations are not present in our formulation, because of the inclusion 
of the "thermodynamic state" equations. 

The obtained complete rate field equations consist of the principles of 
the rate of the virtual power and the rate of conservation of energy, of the 
constitutive relations, and of boundary and initial conditions. These 
formulations provide a sound basis for the formulation of the adopted finite 
element solution procedures. 

One of the most challenging aspects of finite strain formulations is to 
locate analytical solutions with which to compare the proposed formulation. 
Typically, as a first problem, a large strain uniaxial test case was ana- 
lyzed. The case considered examines the rate-dependent plastic response of a 
bar to a deformation history that includes segments of loading, unloading, 
and reloading, each occurring at varying strain and temperature rates. 
Moreover, it was shown that proposed formulation generates no strain energy 
under a pure rigid body rotation. These are surely important demonstrations 

but they only represent a partial test because the principal stretch direc- 
tions remain constant. Finally, a problem which was discussed by Nagtegaal 

and de Jong, and others too, as a problem which demonstrates limitations of 

the consitutive models in many strain formulations, is the Couette flow 
problem. This problem is solved as a third example. The results of these 
test problems show that: 

The formulation can accommodate very large strains and rotations. 

The formulation does not display the oscillatory behavior in the 
streses of the Couette flow problem. 

The model incorporates the simplifications associated with rate- 
insensitive elastic response without losing the ability to model 
rate temperature dependent yield strength and plasticity. 

The problem of buckling of shallow arches under transient thermomechani- 
cal load was investigated next. 

The quasi-1 inear nature of the principle of the rate of virtual power 
suggests the adoption of an incremental approach to numerical integration 
with respect to time. The availability of the field formulation provides 
assurance of the completeness of the incremental equations and allows the use 
of any convenient procedure for spatial integration over the domain V. In 
the present instance, the choice has been made in favor of a simple first 


order expansion in time for the construction of incremental solutions from 
the results of finite element spatial integration of the governing equations. 

The procedure employed permits the rates of the field formulation to be 
interpreted as increments in the numerical solution. This is particularly 
convenient for the construction of incremental boundary condition histories. 

Even under the condition of static external loads and slowly growing 
creep effects, the presence of snap-through buckling makes the inertia 
effects significant. In dynamic analyses, the applied body forces include 
inertia forces. Assuming that the mass of the body considered is preserved, 
the mass matrix can be evaluated prior to the time integration using the 
initial configuration. 

Finite element solution of any boundary-value problem involves the 
solution of the equilibrium equations (global) together with the constitutive 
equations (local). Both sets of equations are solved sinultaneously in a 
step by step manner. The incremental form of the global and local equations 
can be achieved by taking the integration over the incremental time step 
At=tj^j^-tj. The rectangular rule has been applied to execute the resulting 
time integration. 

Clearly, the numerical solution involves iteration. A simplified 
version of the Riks-Wempner constant-arc-length method has been utilized. 
This iteration procedure which is a generalization of the displacement 
control method also allows to trace the nonlinear response beyond bifurcation 
points. In contrast to the conventional Newton-Raphson techniques, the 
iteration of the method takes place in the velocity and load rate space. The 
load step of the first solution in each increment is limited by controlling 
the length ds of the tangent. Either the length is kept constant in each 
step or it is adapted to the characteristics of the solution. In each step 
the triangular-sized stiffness matrix has to be checked for negative diagonal 
terms, indicating that a critical point is reached. 

The analysis was performed with the aid of 24 paral inear isoparametric 
elements. The paral inear isoparametric element is intended for the analysis 
of oriented structures where the geometry is such that the thickness is small 
compared to other dimensions. The characteristics of the element are defined 
by the geometry and interpolation functions, which are linear in the thick- 
ness direction and parabolic in the longitudinal direction. Consequently, 


the element allows for shear strain energy since normals to a mid-surface 
before deformation remain straight, but not necessarily normal to the mid- 
surface after deformation. 

The developed solution scheme is capable of predicting response which 
includes pre- and post-buckling with thermal creep and plastic effects. The 
solution procedure was demonstrated through several examples which include 
both creep and snap-through behavior. 

The last set of problems which are under investigation consists of creep 
or thermal buckling with plastic effects of shells of revolution. 

To develop geometrically nonlinear, doubly-curved finite shell elements 
the basic equations of nonlinear shell theories have to be transferred into 
the finite element model. As these equations are in general written in 
tensor notation, their implementation into the finite element matrix formula- 
tion requires considerable effort. 

The nonlinear element matrices were derived directly from the incremen- 
tally formulated nonlinear shell equations, by using a tensor-oriented 
procedure. A modified version of the classical thin shell theory based on 
the Kirchoff-Love hypotheses, capable of large strains and rotations, is 
presently employed. For this formulation, we are using five "natural" 

degrees of freedom per mid-surface shell node: three incremental velocities 

1 2 

and two rates of rotations about the material coordinates a and a . 

This element was introduced to the solution procedure described earlier. 
Few examples of creep and thermal buckling of a cylindrical panel and spher- 
ical cap are currently under investigation. 

In connection with the progress to date, two papers were published in 
the AIAA Journal and one in the "NASA-U. of Akron Sponsored Meeting" Pro- 
ceedings. Copies of these papers have been mailed to the sponsor. Moreover, 
two papers appeared in the Proceedings of the 28th AIAA/ASME/ASCE/AHS SDM 
Conference. Copies of these two papers are attached herewith. 

FUTURE TASKS 

With regard to additional future tasks, one must include the development 
of a finite element accommodating the more general shell theory formulations 
(A,B, & C), and the incorporation of this into a code. Moreover, three- 
dimensional effects, especially those associated with normal stresses, will 


be incorporated, either through modification of the developed shell theory or 
through local-global finite element procedures in which the local part will 
be based on three-dimensional analysis. 

Finally, another important extension is to replace the material consti- 
tution (initially homogeneous and isotropic-metal) by a more general one that 
allows the analysis of layered fiber-reinforced composites. 
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SOLUTION METHODS FOR ONE-DIMENSIONAL VISCOELASTIC PROBLEMS 

by 

John M. Stubscad* and George J. Simitses** 
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Abstract 

A recently developed differential methodology for solution 
of one-dimensional nonlinear viscoelastic problems is presented. 
Using the example of an eccentrically loaded cantilever beam- 
column, the results from the differential formulation are com- 
pared to results generated using a previously published integrttl 
solution technique. It is shown that the results obtuned from 
these distinct methodologies exhibit a surprising high degree of 
correlation with one another. A discussion of the various {actors 
affecting the numerical accuracy and rate of convergence of these 
two procedures is also included. Finally, the influences of some 
“higher order” effects, such as straining along the centroidal axis 
are discussed. 

I. Introduction 

A number of methods are avmlable to solve viscoelastic 
problems in which the behavior of the material may be ade- 
quately characterized by a linear viscoelastic operedor and where 
the deformation of the body is sufficiently small to allow the 
use of a linear kinematic formulation [1,2]. Commonly, integral 
transform methods, separation of variables, series expansions 
or other such techniques provide methodologies wherein exact 
closed form solutions may be derived. When exact solutions 
cannot be obtained, approximate techniques, such as one pro- 
posed by Schapery [3], may provide an alternate approach. 

The inclusion of nonlinear effects in the analysis signif- 
icantly reduces the mathematical tractability of the problem. 
These nonlinear influences can be induced by geometric factors 
resulting from the magnitude of the deformation or from gross 
rotation of cross sections. Alternatively, nonlinearities in the 
material response may need be included to provide an accurate 
model for material behavior. 
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Independent of whether these nonlinearities are produced 
by geometric or material effects, they invariably result in an 
overall formulation governed by nonlinear equations. Thus, the 
solution methods mentioned above, applicable to linear prob- 
lems, cannot be employed. As such, approximate solution meth- 
ods have been developed and are routinely employed to analyze 
such problems [4j. 

One of these methods is to idetdize the problem is such a 
manner so as to generate inherent simplifications to the govern- 
ing relations. A classical example of this technique was the uti- 
lization of an ideal’ “I” cross sectional geometry in early column 
creep buckling studies [5]. With this approximation, the equa- 
tions governing equilibrium of the column were reduced to much 
simpler forms involving the “average” stresses in the flanges of 
the ideal beam-column. 

Another approach used extensively was to restrict consid- 
eriitions to only certain types of time dependent material be- 
havior [6j. In some cases, this involved retaining only secondary 
creep behavior in the material model. Alternatively, and es- 
pecially when “power law” type constitutve laws were used, 
the constants or exponents of the law were restricted to spe 
cial values for which closed form solution was possible [7j. In 
a few cases, this simplification of the material model, as well 
as the aforementioned geometric simplification technique, were 
employed simultaneously to enable solution. A survey of most 
of these methods has been provided by Hoff [8j. 

An exact numerical solution technique for geometrically 
nonlinear viscoelastic problems has been presented by Rogers 
and Lee [9]. A recent paper by the authors [10] provides a 
method for bounding the solution for such problems. In both 
of these techniques, the solution to the problem is formulated 
in terms of an integral equation which is nonlinear with re- 
spect to both time and space. FVom this, solutions may then 
be readily obtained using relatively standard numerical tech- 
niques. Generally, both the exact and bounding technique can 
be employed for problems wherein the response of the mate- 
rial may be adequately chttracterized using a linear viscoelastic 
model, but where the resulting time dependent deformation of 
the body warrants the use of a nonlinear kinematic formulation. 
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Since the integral solution technique can be applied only 
when the response of the niaterial either is or may be reason- 
ably approximated as being linear, problems involving nonlinear 
viscoelastic material behavior cannot be addressed using these 
methods. Unfortunately, many materials, and especially the ele- 
vated temperature behavior of most metals, require such nonlin- 
ear characterizations. Consequently, an alternate solution proce- 
dure for one-dimensional problems involving nonlinear kinematic 
and nonlinear material effects has been developed. This method, 
hereinafter referred to as the differential formulation, is based 
on the direct solution of the nonlinear differential equations of 
equilibrium. 

Similar to the integral method, the overall procedure is 
predicated on the assumption of a quasistatic response. This 
assumption effectively “decouples” the temporal and spatial de- 
pendence of the problem in such a manner so as to allow the 
general solution to be treated as the combination of the solu- 
tions to a nonlinear “boundary value” problem and a nonlinear 
• initial value" problem. The first of these, the equations char- 
acterizing the time dependent states of quasistatic equilibrium, 
are solved through the use of a Newton type method [11]. In 
contrast, the ’■initial value” problem, resulting from the non- 
linear constitutive law, governs the manner by which the body 
progresses from one state of quasistatic equilibrium to the suc- 
ceeding one. Numerical solutions for this part of the problem are 
generated using a fourth order Runge-Kutta method. Using this 
technique, problems involving the nonlinear thermoviscoelastic 
!)ehavior of thin structural members have recently been exam- 
ined 12. 

In addition to presenting the differential formulation tech- 
nique. a comparison of results obtained through the use of the 
integral and differential formulations is provided. The problem 
of an eccentrically loaded viscoelastic cantilever beam-column is 
employed as the vehicle through which the comparison is per- 
formed. Because of the inherent limitation of the integral tech- 
nique. this comparison is restricted to the consideration of a 
linear viscoelastic material. The specific case considered is that 
of the three parameter viscoelastic solid which has been exam- 
ined in a number of studies [10,13,14]. The results obtained 
from these. two distinct methods of solution exhibit a surpris- 
ingly high degree of correlation with one another thereby estab- 
lishing a high level of confidence in the validity of the methods. 
Finally, the differential formulation is employed to examine the 
influences of some of the “higher order” effects for the class of 
problems under consideration. 

II. Integral Formulation 

Since the integral solution is available in literature [9,10], 
only a brief outline of its development is presented herein. The 


solution is based on the assumption that the beam-column is 
thin and composed of a linearly viscoelastic material. Its ge- 
omerty in the deformed configuration is illustrated in Fig. 1. 
Note that the eccentric load is assumed to be applied quasi- 
statically and its direction does not vary with time. 

Reference line extensional strains are assumed to be neg- 
ligibly smalL Thus the coordinate s' is employed to specify po- 
rtion in both the initial and deformed configurations. A nondi- 
mensional coordinate, s, is defined by dividing s' by the length 
of the beam, L. Assuming a linear distribution of the strains 
through the depth, bending thus occuring within a Bernoulll- 
Euler context, results in a moment-curvature relationship given 
by 




where K{s',t) denotes the curvature and Af(s',T) the bending 
moment at location s'. I is the moment of inertia of the beam 
and J{t) the creep compliance of the material. For the eccentric 
load, R(l), the moment at position s' is given by 


M(s',r) = R[t][S{t) -b acosa(T) - p(s',T)] (2) 



Figure 1. Beam-Column Geometry for Integral Formulation 


From kinematic considerations, it is noted that 


«(s',l) = 
dx(s',t) 


ds' 

dy(s\t) 

ds' 


ds' 

= COS(^{s',f) 


= sin^(s',() 


The boundary conditions for the problem are 


(3) 


^( 0,0 = 0 

(4) 

M{/ -) = aR[t) cos <f>(L,i) 
where a{t) = <p[L,t) for a “rigid” extension. 

Substituting Eqns. (2) and (3) into Eqn. (1) followed by 
differentiation with respect to s yields the governing differen- 
tiid equation. Using the methodology detailed in [9], it can be 
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shown that, for quiescent Initial conditions, the solution for this 
equation is given by 


<>(4,0 = (y ) [,/(O)iZ{/)0(4,/) + J'U - r)fl(r)0(i,T)dT| 


(5) 


where 


5u d(^ 

N = ^ + ^(acos<?(Z) + w(i) - (lla) 


and 


M = F(a cos d(L) 4- w(^) - w(£j) 


where N and M, the force and moment resultants, respectively, 
are defined as 


9(4,f) = (y) cos4>(l,0 + j ff(4,u)sin^(4,/)du 


( 6 ) 


and 


( 7 ) 


g(4,u) = u, 0 < ti < 4 

= 4, 4 < U < 1 

and where the prime, { ' ), in Eqn. (5) denotes differentiation 
with respect to the argument of the function. 

Note that Eqn. (5) represents a time convolution of a non- 
linear Fredholm type integral equation, Eqn. (6). Numerical so- 
lutions for this set of equations are accomplished using Picard’s 
method of successive substitutions [15]. The spatial integrals 
are approximated using Newton-Coates formulae; a fixed step 
trapezoidal rule is used for the time convolution. 


III. Differential Formulation ■ 

Similar to the integral formulation, the differential formu- 
lation for this problem is based on the assumption that bending 
of the beam occurs in accordance with the Euler- Bernoulli hy- 
potheses. The functions u(s,t) and w(4,f) are employed to de- 
note, respectively, the axial and transverse deflection of points 
on the centroidal axis. From this, the extensional strain aJong 
the centroidal axis, c„, is approximately given by 


5u 1,3 w.2 


( 8 ) 


Note that the term 


has been neglected as small in com- 
parison to 1^ . If, in addition, both the strain at the centroidal 
axis and are small in comparison to 1 then it is relatively 
simple to show that 


( 9 ) 


d<t> 9w 3*u 

ds ds ds^ ds- 

Thus, for a linear variation of strain through the cross section, 
this yields 


d<t> 


(10) 


Employing the Principle of Virtual Work followed by inte- 
grating by parts and then manipulating algebraically eventually 
yields the equilibrium equations 


N = y 0)1 d.4 and M = y qon d.4. (12a.(i) 

Based upon an additive decomposition for total strain. Cf 
such that e, = e, -I- c* where <, and represent the elastic and 
creep strain components, respectively, yields, after substitution 
into Eqns. (11) and (12), 

EAco = -f (l + |j + ^(<»tos<>(£)-»’(4)-w(Z))) -Nc (13o) 
and 


El— = F(ocosc>(£,) T w(s) - w(I)) - Mf (13f>l 
ds 

where A represents the area of the cross section. / denotes its 
moment of inertia and the ‘‘pseudo-resultants" Se and are 
defined by 


Nr = ^ EicdA and M, = J r)E(,dA (14a. 6). 

Numerical solutions for Eqns. (14) are computed using a 
modified Newton type method suggested by Thurston 15 . To 
illustrate this method consider a nonlinear differential term of 
the form du^dw" where du tmd dw are differentials of the func- 
tions u and w and m and n represent integer exponents. Assum- 
ing that close trial solutions w and u are available which differ 
from the exact solution by the small quantities Aw and Au so 
that u = u -f Au and w = w Aw, then 


du’" dw" =du’" dw" -h mdu’""' dw" d(Au) 

-I- ndu"* dw""’ d(Aw) -(- /(u. w)Oi Au. Aw 

where /( ) denotes a nonlinear function of u and w and O' 
indicates terms of order AuAw and higher. If the trial solution 
is indeed close to the true solution, then the corrections. Au 
and Aw, will be small. Consequently, the quadratic and higher 
order terms in the corrections will be negli'-'' le in comparison 
to the linear terms and therefore may be neglected. Thus, under 
these conditions, the right-hand side of Eqn. (15) may be closely 
approximated by the linearized form consisting of just the first 
three terms on the right-hand side. 
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With this type of procedure, the original nonlinear differ- 
ential equation is incrementally approximated by a linearized 
form. Employing standard finite-difference formulae, the lin- 
earized form is then converted into a system of tdgebraic equa- 
tions where the unknowns are the corrections to the trial solution 
at the nodes of the finite difference mesh. These relations are 
solved for these corrections, the trial solution is adjusted and 
the process repeated until convergence is obtained. 

.Numerical solution of Eqns. (13) also requires evaluation 
of the “pseudo-resultants” Nc and Me at each point of the finite- 
difference mesh. This is accomplished by evaluating the accu- 
mulated creep strain at select number of points across the cross 
■ section and then using a three perint Newton-Coates quadrature 
formula repetitively to approximate the area integral. Evalua- 
tion of the accumulated creep strain at each of these points is 
accomplished through the use of a fourth order Runge-Kutta in- 
tegration routine to integrate the constitutive law expressed in 
a differential format. 


IV. Example Problem 

The specific example considered is that of a 12 in. long 
beam-column. For simplicity, a square cross section of dimension 
0.5 in. has been assumed. It is also assumed that the beam- 
column is fabricated from a material which can be modeled as 
a three parameter viscoelastic solid. The creep compliance for 
this model, illustrated in Fig. 2, is given by 


.f(<)/J(0) = l-)-!E,/E,Ie-'/'* (16) 

where r, = pi/Ei. For these computations, the numerical values 
for the paramolers have been selected so that = 1. Thus, 
integer values for time equal multiples of the time constant of 
the material. 
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Figure 2. Ideal Three Element “Limited" Creep Model. 


A five point grid in the transverse direction is used for 
computation of the “pseudo-resultants” in the differential for- 
mulation. The points are equidistantly spaced with the first 
and last located at the extreme fibers and the central point po- 
sitioned on the centroidal axis. 

Since the solutions of both the' differential and integral for- 
mulations are only satisfied at a discrete number of points over 
the length of the beam-coltunn, the first question that must be 


addressed is how sensitive the results are to the number of points 
used in the approximation. Table 1 provides a comparison of the 
variation of the results from the differential formulation for the 
initial elastic deflection of the beam- column as the number of ap- 
proximating points is doubled from 10 to 20 and then doubled 
again to 40 for an eccentricity ratio of 0.05 and an applied load 


Table 1. Differential Elastic Solution Versus Number 
of Nodes for P/P, = 0.75 and a/L = 0.05 


TVcinsverse Deflection for Various Numbers of Elements (in) 


Number of Elements 


5 = s'/L 

40 

20 

% Diff 

0.0 

0.000000 

0.000000 


0.1 

0.024673 

0.024705 

0.130 

0.2 

0.098189 

0.098307 

0.120 

0.3 

0.219074 

0.219332 

0.118 

0.4 

0.384931 

0.385376 

0.116 

0.5 

0.592512 

0.593187 

0.114 

0.6 

0.837815 

0.838755 

0.112 

0.7 

1.116191 

1,117426 

0.111 

0.8 

1.422456 

1.424007 

0.109 

0.9 

1.750997 

1.752879 

0.107 

1.0 

2.095877 

2.098098 

0.106 

s = s'/L 

40 

10 

% Diff 

0.0 

0.000000 

0.000000 


0.1 

0.024673 

0.024809 

0.551 

0.2 

0.098189 

0.098763 

0.585 

0.3 

0.219074 

0.220289 

0.555 

0.4 

0.384931 

0.387110 

0.566 

0.5 

0.592512 

0.595755 

0.547 

0.6 

0.837815 

0.842435 

0.551 

0.7 

1.116191 

1.122177 

0.536 

0.8 

1.422456 

1.430090 

0.537 

0.9 

1.750997 

1.760165 

0.524 

1.0 

2.095877 

2.106817 

0.522 


{% differences are with respect to 40 clement solution) 


of P/P, = 0.75. The Euler load, P,, here is based on a perfect 
geometry and use of the instantaneous compliance of the mate- 
rial, J(0). Table 2 provides a similar comparison for the szime 
loading and eccentricity ratio but for the integral formulation. 

From these two Tables it can be seen that there is very lit- 
tle change in the computed transverse deflection as the number 
of approximating points is increased. In both cases, the initial 
elastic solution of the ten element model is within 1.0 % of the 40 
element model results. Additionally, the relative magnitude of 
the errors between the 10 and 40 element models of the differ- 
ential formulation are virtually identical to those exhibited by 
the integral solution methodology. Although not specifically in- 
dicated on these Tables, it should also be noted that the angles 
of rotation are not always smtdL For example, at this load and 
eccentricity, the initial elastic end rotation of the beam-column 
exceeds 17 deg. Of course, smaller rotations are exhibited at the 
lower loads and lower eccentricities. 
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Table 2. Integral Elastic Solution Versus Number 
of Nodes for P/P, = 0.75 and ajL = 0.05 


Transverse Deflection for Various Numbers of Elements (in) 


Number of Elements 


s = s'/i 
0.0 

40 

0.000000 

20 

0.000000 

% Diff 

0.1 

0.024616 

0.024652 

0.146 

0.2 

0.097967 

0.098125 

0.161 

0.3 

0.218586 

0.218890 

0.139 

0.4 

0.384079 

0.384691 

0.159 

0.5 

0.591200 

0.592089 

0.150 

0.6 

0.835951 

0.837267 

0.157 

0.7 

1.113686 

1.115374 

0.152 

0.8 

1.419227 

1.421422 

0.155 

0.9 

1.746972 

1.749594 

0.150 

1.0 

2.091003 

2.094171 

0.152 

s = s'/i 
0.0 

40 

0.000000 

10 

0.000000 

% Diff 

0.1 

0.024616 

0.024781 

0.670 

0.2 

0.097967 

0.098753 

0.802 

0.3 

0.218586 

0.220123 

0.703 

0.4 

0.384079 

0.386547 

0.643 

0.5 

0.591200 

0.595443 

0.718 

0.6 

0.835951 

0.841727 

0.691 

0.7 

1.113686 

1.120936 

0.651 

0.8 

1.419227 

1.429034 

0.691 

0.9 

1.746972 

1.758422 

0.655 

1.0 

2.091003 

2.105685 

0.702 

(% differences 

are with respect to 40 element 

solution) 


Therefore, it is concluded that both formulations exhibit 
the same low level of sensitivity to the number of elements used 
in the analysis. Concurrently, these results also indicate that 
a 10 element model crin be used with either method without 
generating significant errors in the analysis. It should be noted 
that all of the differential formulation results presented above 
are based on the use of an exact expression for evaluating the 
angle of rotation. 

A direct comparison between the results generated by the 
two formulations is provided in Tables 3 and 4. Again, the com- 
parison is based on the 0.50 eccentricity ratio which inherently 
produces larger angles of rotation. The angle of rotation of the 
cross section, determined from the integral technique, are also 
provided. 

Table 3 is based on the use of an exact expression for 
evaluation of sin^ in the solution of the differential methodol- 
ogy. The “approximate” sin^ results, provided in Table 4, on 
the other hand, are based on use of the common approxima- 
tion sin<^ a: -dvr/ds, to calculate the angle of rotation. Note 
that, other than this particular change, these two differential 
formulations are, otherwise, completely identical. Both include 
the influences of axial as well as transverse deflection on the be- 
havior of the beam-column. The influence of deleting the axial 
deformation on the accuracy of the predicted results is discussed 
later. 


Table 3. Integral and Differential ( ^ ) Elastic 

Solutions for Various Loads aJL = 0.05 


Transverse Deflection from Various Solutions 

(in) 



Exact 





sin 6 


Angle 

s = s' JL 

Integral 

didt 7c 

Diff 

(deg) 


for P/P, = 0.25 



0.0 

0.000000 

0.000000 


0.00 

0.1 

0.002615 

0.002616 

0.04 

0.25 

0.2 

0.010450 

0.010448 

-0.02 

0.50 

0.3 

0.023444 

0.023445 

0.00 

0.74 

0.4 

0.041520 

0.041531 

0.03 

0.98 

0.5 

0.064590 

0.064589 

0.00 

J oo 

0.6 

0.092475 

0.092483 

0.01 

1.44 

0.7 

0.125010 

0.125031 

0.02 

1.66 

0.8 

0.162039 

0.162044 

0.00 

1.87 

0.9 

0.203256 

0.203279 

0.01 

2.07 

1.0 

0.248513 

0.248497 

-O.OI 

2.25 


for P/P, = 0.50 



0.0 

0.000000 

0.000000 


0.00 

0.1 

0.008283 

0.008285 

0.02 

0.79 

0.2 

0.033056 

0.033044 

-0.04 

1.57 

0.3 

0.073940 

0.073948 

0.01 

2.33 

0.4 

0.130441 

0.130516 

0.06 

3.07 

0.5 

0.202009 

0.201999 

0.00 

3.76 

0.6 

0.287513 

0.287559 

0.02 

4.40 

0.7 

0.385945 

0.386070 

0.03 

5.01 

0.8 

0.496357 

0.496382 

0.01 

0..52 

0.9 

0.616919 

0.617049 

0.02 

6.01 

1.0 

0.746772 

0.746670 

-0.01 

6.38 



for P/P, = 0.75 



0.0 

0.000000 

0.000000 


0.00 

0.1 

0.024781 

0.024809 

0.11 

2.37 

0.2 

0.098753 

0.098763 

0.01 

4.68 

0.3 

0.220123 

0.220289 

0.08 

6.91 

0.4 

0.386547 

0.387110 

0.15 

9.03 

0.5 

0.595443 

0.595755 

0.05 

10.97 

0.6 

0.841727 

0.842435 

0.08 

12.68 

0.7 

1.120936 

1.122177 

0.11 

14.23 

0.8 

1.429034 

1.430090 

0.07 

1.5.41 

0.9 

1.758422 

1.760165 

0.10 

16.44 

1.0 

2.105685 

2.106817 

0.05 

17.04 


{% differences are with respect to integral solution) 
(comparisons based on results from 10 elemeiii iiioileM 


As the data of these Tables indicate, there are no sisTnifi- 
cant differences in predicted transverse deflections. Even for an 
end rotation angle of 17 deg, the differences between the various 
results is well below that required in engineering computations. 
It should be noted that this high level of correlation continues 
to exist for even greater angles of rotation. The reason why 
this high correlation exists is that, under compressive loading, 
the derivative of the axial displacement is negative. Since, the 
square of the slope of the transverse deflection is always positive, 
these algebraically combine so that the sum is lower in magni- 
tude than either of the individual terms. This serves to reduce 
the magnitude of the centroidal axis strain. Because the differ- 
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Table 4 . Integral and Differential {■§;) Elastic 

Solutions for Various Loads aJL = 0.05 


Transverse Deflection from Various Solutions (in) 
Approx. 

sin 6 Angle 

5 = .«';i Integral djdt % Diff (deg) 



for P/P, = 0.25 



0.0 

0.000000 

0.000000 


0.00 

0.1 

0.002615 

0.002616 

0.00 

0.25 

0.2 

0.010450 

0.010448 

-0.02 

0.50 

0.3 

0.023444 

0.023445 

0.00 

0.74 

0.4 

0.041520 

0.041531 

0.03 

0.98 

0.5 

0.064590 

0.064589 

0.00 

1.22 

0.6 

0.092475 

0.092483 

0.01 

1.44 

0.7 

0.125010 

0.125031 

0.02 

1.66 

o.s 

0.162039 

0.162044 

0.00 

1.87 

0.9 

0.203256 

0.203279 

0.01 

2.07 

1.0 

0.248513 

0.248497 

-0.01 

2.25 



for P/P, = 0.50 



0.0 

0.000000 

0.000000 


0.00 

0.1 

0.008283 

0.00S286 

0.04 

0.79 

0.2 

0.033056 

0.033044 

-0.04 

1.57 

0.3 

0.073940 

0.073949 

0.01 

2.33 

0.1 

0.130441 

0.130516 

0.06 

3.07 

0.5 

0.202009 

0.202000 

0.00 

3.76 

0.6 

0.287513 

0.287560 

0.02 

4.40 

0.7 

0.385945 

0.386071 

0.03 

5.01 

0.!> 

0.496357 

0.496383 

0.01 

5.52 

0.9 

0.616919 

0.617050 

0.02 

6.01 

1.0 

0.746772 

0.746672 

-0.01 

6.38 



for P/P, = 0.75 



0.0 

0.000000 

0.000000 


0.00 

0.1 

0.024781 

0.024310 

0.12 

2.37 

0.2 

0.098753 

0.098767 

0.01 

4.68 

0.3 

0.220123 

0.220299 

0.08 

6.91 

0.4 

0.386547 

0.387128 

0.15 

9.03 

0.5 

0.595443 

0.595782 

0.06 

10.97 

O.C 

0.841727 

0.842473 

0.09 

12.68 

0.7 

1.120936 

1.122228 

0.12 

14.23 

O.S 

1.429034 

1.430155 

0.08 

15.41 

0.9 

1.758422 

1.760244 

0.10 

16.44 

1.0 

2.105685 

2.106911 

0.06 

17.04 


I'" differences tire with respect to integral solution) 
(comparisons based on results from 10 element^ models) 


ence between the exact and approximate expressions for sin^ is 
related to the y'l + 2«<, in the denominator of the exact expres- 
sion. reducing the magnitude of the centroidal strain inherently 
improves the accuracy of the approximation. 

This influence is best illustrated by the data of Table 5. 
Here, the integral solution is compared to an approximate solu- 
tion for which the effect of the centroidal axis strain terms has 
been suppressed. This was accomplished by first eliminating 
a.' the {dvr/ds)^ terras from the governing equations. The EA 
modulus-area product was then artiflcially increased through 
multiplication by a factor of 1000. This second change reduces 
the magnitude of the axial deflection by approximately the same 
factor. 


The net result of these changes is to create a differen- 
tial model where the centroidal axis effectively is inextensional. 
While it might be anticipated that this would improve the cor* 
relation between the differenticd and integral results, such is not 
the case. When the angle of rotation is very smaU, as those 
which result from a low level of loading and a small eccentricity, 
all the formulations provide virtually identical predictions. How- 

Table 5. Integral and Modified Differential (^) 

Elastic Solutions for Various Loads a/I = 0.05 


Transverse Deflection from Various Solutions (in) 


s = s'/L 

Integral 

Modified 
(w/o e,) 
diet 

% Diff 

Angle 

(deg) 

0.0 

0.000000 

for P/P, = 0.25 
0.000000 

0.00 

0.1 

0.002615 

0.002616 

0.00 

0.25 

0.2 

0.010450 

0.010446 

-0.04 

0.50 

0.3 

0.023444 

0.023444 

0.00 

0.74 

0.4 

0.041520 

0.041527 

0.02 

0.98 

0.5 

0.064590 

0.064586 

-O.Ol 

1.22 

0.6 

0.092475 

0.092478 

0.00 

1.44 

0.7 

0.125010 

0.125030 

0.02 

1.66 

0.8 

0.162039 

0.162043 

0.00 

1.87 

0.9 

0.203256 

0.203287 

0.02 

2.07 

1.0 

0.248513 

0.248508 

0.00 

2.25 

0.0 

0.000000 

for P/P, = 0.50 
0.000000 

0.00 

0.1 

0.008283 

0.008294 

0.13 

0.79 

0.2 

0.033056 

0.033073 

0.05 

1.57 

0.3 

0.073940 

0.074032 

0.12 

2.33 

0.4 

0.130441 

0.130665 

0.17 

3.07 

O.S 

0.202009 

0.202274 

0.13 

3.76 

0.6 

0.287513 

0.287975 

0.16 

4.40 

0.7 

0.385945 

0.386711 

0.20 

5.01 

0.8 

0.496357 

0.497264 

0.18 

5.52 

0.9 

0.616919 

0.618270 

0.22 

6.01 

1.0 

0.746772 

0.748235 

0.20 

6.38 

0.0 

for P/P, = 0.75 
0.000000 O.OOOGOO 


0.00 

0.1 

0.024781 

0.025543 

3.07 

2.37 

0.2 

0.098753 

0.101698 

2.98 

4.68 

0.3 

0.220123 

0.227057 

3.15 

6.91 

0.4 

0.386547 

0.399299 

3.30 

9.03 

0.5 

0.595443 

0.615237 

3.32 

10.97 

0.6 

0.841727 

0.870876 

3.46 

12.68 

0.7 

1.120936 

1.161483 

3.62 

14.23 

0.8 

1.429034 

1.481683 

3.68 

15.41 

0.9 

1.758422 

1.825548 

3.82 

16.44 

1.0 

2.105685 

2.186716 

3.85 

17.04 


(% differences are with respect to integral solution) 
(comparisons based on results from 10 element models) 


ever, as the angle of rotation increases, as those which occur at 
the higher loadings and the greater load eccentricity, the mod- 
ified differential formulation predictions begin to diverge from 
those of the other two. The difference between the predicted 
results is most notable for the highest load and largest load ec- 
centricity example. 
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This behavior results from the way the end of the beam- 
column deflects. In the modified model, the end of the beam- 
column basically moves only in the vertical direction (see Fig. 1). 
The standard model, however, allows the end to move both ver- 
tically and horizontally. Thus, for a given vertical deflection, the 
leftward movement predicted by the standard model increases 
the angle of rotation. This reduces the moment by reducing the 
moment arm. Thus, to support a given load, the standard for- 
mulation beam-column model need not deflect as much as the 
modified model must. 

Comparisons of the viscoelastic predictions provided by 
the integral and the exact differential are provided in Table 6. 
A comparison with the approximate sin^ differential solution 
is not included because the differences between the exact and 
approximate results ag^ are so small as to be negligible. These 
computations are based on the use of an ideal three element 
“limited” material model illustrated in Fig. 2. 

As demonstrated by the data of this Table, the high cor- 
relation between the elastic solutions provided by the integrrd 
and differential solution methods carries over directly into the 
viscoelastic analysis. 

The final item meriting discussion is the length of the time 
increment used in each of the formulations. Unlike the prior re- 
sults, some differences do exist between the maximum allowable 
time step increments for the integral and differential formula- 
tions. Additionally, the allowable time step increment for the 
integral formulation exhibits a higher dependence on the actual 
angle of rotation than does the differential formulation. 

In general, a relatively small time step increment must be 
used with the integral solution methodology. For example, the 
results presented above are based (n\ the use of a 0.01 time step. 
As the length of this time step is increased, the accuracy of the 
solution decreases and tends to underpredict the deflection. This 
convergence “from below” is not surprising since the convolution 
integral is approximated as the sum of a finite number of terms. 

In contrast, much larger time step increments can be used 
with the differential formulation. This is principally attributed 
to the high accuracy provided by the Runge-Kutta integration 
routine. Most of the results provided above were developed using 
a 0.10 time step. Additionally, the allowable length for this time 
step increment tended to be rather insensitive to the angle of 
rotation. The allowable time step for the integral formulation, 
on the other hand, exhibited a high level of sensitivity to the 
angle of rotation. Larger angles of rotation required significantly 
shorter time steps for accurate results to be obtained. 

These factors combine in a rather interesting manner with 
regard to which method of analysis is computationally more ef- 
ficient. Typically, for the analysis of short periods of viscoelas- 
tic deformation the integral solution method was two to three 
times faster than the differential method. This is attributed to 


Table 6. Integral and Differential Viscoelastic 
Solutions for Three Time Constants 
P/P, = 0.50, and a/I = 0.05 



Transverse Deflection 

(in) 




Exact 





sind 


Angle 

II 

Integral 

djdt % Diff 

(deg) 


at 

time = 0.0 



0.0 

0.000000 

0.000000 


0.00 

0.1 

0.008283 

0.008285 

0.02 

0.79 

0.2 

0.033056 

0.033044 

-0.04 

1.57 

0.3 

0.073940 

0.073948 

0.01 

2.33 

0.4 

0.130441 

0.130516 

0.06 

3.07 

0.5 

0.202009 

0.201999 

0.02 

3.76 

0.6 

0.287513 

0.287559 

0.02 

4.40 

0.7 

0.385945 

0.386070 

0.03 

5.01 

0.8 

0.496357 

0.496382 

0.01 

5.52 

0.9 

0.616919 

0.617049 

0.02 

6.01 

1.0 

0.746772 

0.746670 

-0.01 

6.38 


at 

time = 1.0 



0.0 

0.000000 

0.000000 


0.00 

0.1 

0.015083 

0.015089 

0.04 

1.44 

0.2 

0.060151 

0.060128 

-0.04 

2.86 

0.3 

0.134319 

0.134344 

0.02 

4 22 

0.4 

0.236430 

0.236612 

0.08 

5.. 54 

0.5 

0.365200 

0.365193 

0.00 

6.76 

0.6 

0.518035 

0.518161 

0.02 

7.86 

0.7 

0.692C34 

0.692946 

0.05 

8.87 

0.8 

0.886869 

0.886967 

0.01 

9.70 

0.9 

1.096666 

1.097006 

0.03 

10.44 

1.0 

1.320154 

1.320004 

-0.01 

10.95 


at time = 2.0 



0.0 

0.000000 

0.000000 


0.00 

0.1 

0.019127 

0.019138 

O.OC 

1.83 

0.2 

0.076254 

0.076232 

-0.03 

3.62 

0.3 

0.170149 

0.170202 

0.03 

5.34 

0.4 

0.299200 

0.299484 

0.09 

7.00 

0.5 

0.461621 

0.461664 

0.01 

8.52 

0.6 

0.653844 

0.654093 

0.04 

9.88 

0.7 

0.872716 

0.873246 

0.06 

11.13 

0.8 

1.115343 

1.115622 

0.03 

12.12 

0.9 

1.376227 

1.376871 

0.05 

12.99 

1.0 

1.652861 

1.652902 

0.00 

1 3. .55 


(% differences are with respect to integral solution) 
(comparisons based on results from 10 element models! 


two factors. The first is the comparitively slow Runge-Kutta 
integration procedure used in the differential formulation. For a 
short period of viscoelastic deformation, the calculation of the 
convolution integral of the integral formulation, requiring sim- 
ple summation of a limited number of terms, can be performed 
much more rapidly. 

The second factor is that the fixed point iteration scheme 
of the integral formulation, although requiring more iterations 
than the Newton method, is also performed more rapidly since it 
is simply an algebraic operation. The Newton method, in con- 
trast, requires inversion of the matrix premultiplying the vec- 
tor of trial function corrections and then numerical evaluation 
through solution of the system of equations. Even for just a 10 
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.■ipiiienr h^ani. this process is slow in comparison to the fixed 
point il<'ralion. 

However. a> the length of the period of viscoelastic defor- 
mation increases, this relative speed relationship reverses. Even- 
tually. the differential formulation begins to generate solutions 
more ranidlv than the integral method. In the example problem 
(iescribeil above, this generally occurred approximately between 
the second and third time constants. The reason for this change 
is directlv related to the computation of the convolution inte- 
cral of the integral technique. As time increases, the number 
of terms in the summation increases linearly. This in turn, in- 
creases the number of additions which must be performed and 
therefore linearly increases the time need for each computation. 
In contrast, the speed of the Runge-Kutta integration routine is 
virtually independent of time. Thus, the continually increasing 
computation effort required in the integral technique eventually 
e.xces-d' that need for the differential technique thereby reversing 
the r- laiive speed relationship. 

V. Conclusions 

Ba-ed on the results reported herein and elsewhere [12] 
.t i- concluded that the differential formulation procedure pre- 
•ented can be employed for the analysis of quasistatic nonlinear 
one-dimensional viscoelastic problems. This conclusion is based 
directiv on the high level of correlation between results devel- 
oped u-ins. this formulation technique to those obtained with 
the previously published integral method for solution of such 
probjem.-. .Additionally, it is noted that both of these meth- 
od- exi.ibii exceptionally similar accuracy characteristics with 
regard to the numher of elements employed in the approxima- 
tion. For both, a relatively low number of elements can be used 
without engendering any significant errors. 
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Abstract 

The proMeiii of buckling of shalluw arches under transient 
tiiermoiiiocliauical loads is investigated. The analysis is based on 
nonlinear geometric and constitutive relations, and is expressed 
in a rate form. Tlie material constitutive equations are capable 
of rejirexlucing all non-isotheriiial. elasto-viscoplastic character- 
istics. The solution scheme is capable of predicting response 
wliicli iticliides pre and postbuckling with creep and plastic ef- 
fects. The solution procedure is demonstrated through several 
examples which include both creep and snap-through behavior. 


L It itroditctio n 

Elastic sna])-through of low arches under quasi-static loads 
has been the subject of several investigations. The signiRcanre of 
the problem, in so far as it illustrates certain important features 
in tnore cotu]iIicated buckling problems of plates and shells, was 
pointed out by Marguerre (1). who constructed a simplified me- 
chanical model to demonstrate these features. Timoshenko |2) 
obtained an approximate solution to the problem of a low arch 
under a uniformly distributed load. Biereno |3) considered the 
problem of a low parabolic arch loaded laterally at the mid-point 
with a concentrated load. He also investigated snap-through 
buckling of a shallow spherical cap, pinned along its circular 
boitndarv. utider the action of a concentrated load applied along 
the axis of rotational syiiunetry. He presented his approximate 
'ohnions by means of load -deflect ion curves and equations from 
which the critical load could be calculated. 

In ly.'ij. Eung and Kaplan [4] investigated the prohlem of 
low pinned arches of various initial shapes and spatial distribu- 
tions of tlie lateral load. Their results show that a very shallow 
arch snaps-through syiimietricaJly, wliereas a higher arch buck- 
les asyituiietrically. They also ran a limited number of experi- 
memal tests, and their experimental data are in good agreement 
with ilirir tiieoretic.il results. About the same time. Hoff and 
Bruce .i'. in investigating the possiliility of snap-through buck- 
litig of a low pinned arch witli a half-sine-wave initial shape 
ttnder a half-sine-wave distributed dynamic load (suddenly ap- 
plied with iiiHtiite duration), obtained results for the quasi-static 
load case which are identical to tliose obtained by Fung and Ka- 
plan for file same problem. 
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In 1902, Gjelsvik and Bodner [C] obtained an approximate 
solution to the problem of a low circular arcli with a concen- 
trated load at the mid point of tlie arch and clamped bound- 
ary conditions. They also reported on experimentai tests, and 
there is good agreement between their theoretical jjredictioiis 
and their exjtprimental results. Schreyer and Masur |7i obtained 
an exact solution to their problem (and for the load case of 
uniform pressure), and they showed that for the concentrated 
load case.the arch snaps-through syimnetricallv regardless of 
tlie value of the rise parameter. Masiir <ind Lo '8 presented 
a genera] discussion of the behavior of the shallow circular arch 
regarding buckling, post -buckling and imperfection sensiiivitv. 
Snapping of tow pinned arches resting on an elastic roundation 
has been investigaletl by .Simitses [9[. This system exhibits ,il| 
forms of oxi>erimentaUy observed buckling phenomena (smooth 
and violent ) and of theoretically predicted responses (limit point, 
bifurcation with stable branching and bifurcation with unstable 
branching), and it is presented with sufficient detail in Ref. 10. 
Experimental results have also been reported by Roorda ;11]. 

The effects of inelastic material behavior found tlteir wav into 
the literature since the 19C0*s. Onat and Shu jl2; considered the 
behavior to l>e one of rigid-perfectly plastic. Fromciosi. August! 
and Sparacio (13) discussed the collapse of arches under repeated 
loads with inelastic iimterial behavior. Studies of inelastic snap- 
through buckling of shallow arches were also reported by Lee 
and Murphy |14). In addition August! [15] investigated plastic 
buckling of a model of a three hinged arch in 19G8. and a tnore 
complete analysis of the same model was provided by Batter- 
man jlOj in 1971. Finally, the reader who is interested in the 
ultimate strength of p.ir.ibolic steel arches with bracing system 
is referred to Komatsu |17|. who considers inelastic in-planc and 
out-of-plane instabilities and proviiles ilesign formula ftir each 
rase. 

The elastic response of arches under sudden (dynamic) ap- 
plication of the extental loads has been reported by Hoff and 
Bruce [5). Hsu [21]. l22j. and Lock [23j. For a more complete 
l>it)L'ogra]>liy see Ref. 24. 

Creep buckling of shallow arches has been investigated bv 
Huang and Nachbar [18], Krajeinovic [19]. and Botros and Bi- 
enek |20]. 

As far as the authors know no work ha< been reported on the 
non-isothermai elastovisro-plaslic behavior of shallow arches. 
The purpose of this paper is to demonstrate the effect of Iiighly 
nonlinear material behavior on the snap througii and creep Itucfc- 
ling of shallow arches. 


iiS 



11. Eiaslo-Theriiio-Visfoplastic Constitutive Relations 

The prediction of buckling loads and postbuckling behavior 
of structural coiupononts. like shallow arches* made of a realistic 
iiietalljc material and subjected to non-isothermal theniiome- 
chanical loads has increased in importance in recent years. 

ruder this kind of severe loading conditions, the structural 
behavior i' highly nonlinear due to the combined action of ge- 
ometrical and physical nonlinearities. On one side, finite defor- 
mation in a stressed structure introduces nonlinear geometric 
effects. Or, the other side, physical nonlinearities arise even in 
-mal! strain regimes, whereby inelastic phenomena play a par- 
ticulariv important role. From a theoretical standpoint, nonlin- 
ear constitutive equations should be applied only in connection 
with nonlinear transformation measures (implying both defor- 
mation and rotations). However, in almost all of the works in 
this area 2oh tlie two identified sources of nonlinearities are 
always separated. The separation yields, at one end of the spec- 
trum. pro!, loins of large response, while at the other end, prob- 
1-iiis of viscous and. or non-isothermal behavior in the presence 
, ,f small st r.'iiii- 

The cla-'ical theories, in which the material response is char- 
artori.’s-d a combination of distinct elastic, theniial, tiineinde- 
pondoiit iiiolastic (plastic) and time dependent inelastic (creep) 
(lof.rmaiion components cannot explain some phenomena, that 
can lio oti-orved in complex therino-ntechanical loading histo- 
ries. This is particularly true when high-temperature non- 
i'Othertual processes must he taken into account. There is a 
sizeable body of literature (25J, i26l on phenomenological consti- 
tutive equations for the rate and temperature dependent plastic 
deformat ioti behavior of metallic materials. However, almost all 
of these new "unified" theories are based on small strain theo- 
ries Afi,i s.-n-ral suffer from .some thennodvnanuc inconsistencies. 

In a previous paper [27), the authors have presented a com- 
plete of Constitutive relations for nonisothermal. large strain, 
elasio visfoidastic behavior of metals. It was shown there [27] 
that the metric tensor in the converted (material) coordinate 
sv'ieni ran be linearly decomposed into elastic and (visco) plas- 
tic parts. So a yield function was assumed, which is dependent 
on the -ate of change of stress on the metric, on the temperature 
and on a set of internal variables. Moreover, a hypoelastic law 
vsa- rliosen to describe the thermo-elastic part of the defomta- 
iion. 

time and temperature dependent viscoplasticity model was 
formulated in this cons-ected material system to account for finite 
strains and rotations. The history and temperature dependence 
were incorporated through the introduction of infernal variables. 
The choice of these variables, as wed as their evolution, was mo- 
tivated by thermodynamic considerations. 

The nonisothermal elasto-viscoplastic defonuation process 
was described completely by ’’thermodynamic state” equations. 
.Most investigators [25], [26] (in the area of viscoplasticity) ein- 
plov plastic strains as state variables. The authors’ previous 
study 27 shows that, in general, use of plastic strains as state 
variables may lead to inconsistencies with regard to thennody- 
namir considerations. Furthermore, the approach and formu- 
lation employed in previous works leads to the condition that 
all the plastic work is completely dissipated. This, however, 
is ir, contradiction with experimental evidence, from which it 


emerges that part of the plastic work is used lor jiroaucmg l esiu- 
ual stresses in the lattice, which, when phenomenologically con- 
sidered, causes hardening. Both limitations were excluded from 
this |27] formulation. Accuracy of the formulation was checked 
on a wide range of examples [28]. 

The constitutive relation will be rephrased here in some dif- 
ferent form. For brevity we compile only some notations and 
fundamental relations which are used in the formulation of the 
intended constitutive low. For details, see [27] and [28J. 

Conrerning the formulation of constitutive laws it is advan- 
tageous to use a material (co-moving) coordinate system. The 
transformation from the underforined state (metric go,) to the 
deformed state can be represented by the tensor: 

fl=9"9rk or {/‘Ml = S" ffri, (1) • 

The total deformation rate is defined by 

4 = \9’'9.k = -\9ir9'" = lu-'Um = ( 2 ) 

here (') denotes time material derivative. The expression 

Ji= uYk + 4 /; - difi = w 

represents the synunetric part of (/)'^ or the rovariant deriva- 
tive with re.spert to time, also railed the ronverlcd derivation, 
which is due to Zareiuba and Jaumaiin. 

The total deformation ran be decomposed according to 

fl =g"”9,„r g.k =nii (4) 

• 

where the superscript ( ) relates to a fictitious configuration de- 
fined by a fictitious reversi)>le process with frozen internal vari- 
ables. The decomposition of Eq. (4) leads to an additive de- 
composition of the deformation rate 

l») (0 

d\ =di -f 4 (5) 

('} <’.l 

d* is related to the reversible deformations, whereas dj, denotes 
the remaining part of the deformation rate. 

For the description of the stress state, we introduce the Kir- 
chofT stress tensor sj,, which is connected with the real Cauchy 
stress tensor trj,, by the relation; 

O 

4 = ^4 (6) 

Assuming that the elastic behavior is not affected essentially 
by inelastic defonuations, the following hypoelastic incremental 
law was chosen [27] 

4=2jj't+{^4 + «m (7) 

with 

t\ : weighted stress deviator 

G : shear modulus 

K : bulk modulus 

a : coefficient of thermal ex)>aasion 

( 

The following constitutive relations were established [27] for 
the inelastic behavior, 
yield condition: 

F = [ti-cP j4)(f* - c P s4) - h^[A,T] = f^-k^>0 (8) 
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arcoinpaiiving equilibrium stale: 


F = {i[- CP -fpgtii)- I^'HA,T) = /’ - i-* = 0 (9) 


evolution law for inelastic deforniations: 




vith 


and 


d'^.= 2X01 - c P gl3l) 


_ 1 ..I u',-^P3A)U^ _ 

~ ■ 'y 


A 7) 


= 


1 


1 + 4t;A 




xUi- ef>gi3i) + cpgdi 


1) 


evolution laws for the internal variables: 

A = l/I. rff 

P 

r (I I 

if 

dF ^ 


then 


with 


then 


r = 0 nnd 

(O 

4=rfl 

bl ^ .. . 0 . 

0 and </;.= 2A(ll. - c P p,'3L) 

f < 0 


('I 

rfj, = dlt 

/i = 0 

,]■ = 0 


( 10 ) 

(11) 

(12) 

(13) 

(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

( 20 ) 

( 21 ) 


VVitliin iJie developed frame the elaslo-viscojdastic behav- 
ior is governed by the scalar material functions c[sl.T. A,0[). 
l-( .4. r ). g{ sj^, .A,T.4).i(A.T,l3l),!oid 7i{A,T, l5 ], ). These ma- 
terial functions can be determined from a series of monotonic 
and cyclic jirocesses with proportional and nonproportional paths 
at difr''renl temperature levels '(28). 


111. General Formulation and Solution Schemes 

The rale form of the constitutive equations suggests that a 
rate approach lie taken toward the entire problem so that flow 
is viewed as history dependent process rather than an event. 
For this purpose, a complete true ab-initio rate theory of kine- 
matics and kinetics for continuum and curved thin structures, 
without any restriction on the magnitude of the transformation 
was presented in Ref. 28. It is implemented with respect to 


a body- fixed system of converted coordinates, and it is valid 
for finite strains and finite rotations. The time dependence and 
large strain behavior are incorporated through the introduction 
of the time rates of change of the metric and of the curvature. 

The constitutive law may be applied to the conservation of 
momentum via an appropriate variational principle. The prin- 
ciple of virtual power (or of virtual velocities) reads 

^ (r’^6vj,id\’ - pfHvjdV - vVh-jd.A = 0 (22 ) 

where 6vj are the virtual velocities. /■' the bodv forces per unit 
mass and vT^ the surface tractions. Total differentiation of Eq. 
(22) yields, 



At any instant. Eq. (23) must be satisfied. The virtual ve- 
locity and its time derivative are then independent. Moreover, 
the last three terms of Eq. (23) are equivalent to Eip i22). 
Hence, the principle of the rale of virtual power may he ol>- 
tained in its concise form. For further classifications, tlie total 
derivative of the stress components will be represented by the 
Jauman derivative, and the following integrals are defined by 

I, = a’J 6vj,,d\- (241 

Ij = (2.M 

I, = j^^■>yjh■J,,dV (2(ii 

Then, .substitution in Eq. (23) yields the final form of the itrin- 
ciple of the rale of virtual power. 

I = I, + Ij + I, = ~ (27 ) 

The qiiasi-linear nature of the principle of the rate of virtual 
power suggests the adoption of an incremental approach to nu- 
merical integration with respect to time. The availability of the 
field formulation provides assurance of the complelene's of the 
incremental equations and allows the use of any convetiieni pro- 
cedure for spatial integration over the dottiain V. In the present 
instance the choice has been made in favor of a simple lir't order 
expansion in time for the construction of incremental solutions 
from the results of finite element spatial integration of the gov- 
erning equations. 

The procedure eni|)loyed permits the rales of the field formu- 
lation to be interpreted as increments. in the numerical solution. 
This is particularly convenient for the construction of incretiien- 
ta) boundary condition histories. 

The finite element method for spatial discretization has been 
well documented (see. e.g. Zienkiewicz [29; or Oden 30 ) and 
will not be detailed here. Restricting attention to a single finite 
element B„, the velocity field is approximated bv r,(j-'|. 

f,(z"’)'-f-.(z"') = (r“'^)^v'r(j-')l‘* <1.-3 =1 .V (28) 
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hi I'-ij. r ’ arc generalized nodal velocities. F”'’ is depen- 

d-iiT iKKin ilie present nodal coordinates. is a vector of 

[ir'-rr i.er! functions of (z-'). and .V is the number of degrees of 
freedom associated with the elenienl. The matrix F"*^ is defined 
I'v r' puiriiii’ tliat evaluation of r, at nodal positions yields . 

Tile algebraic counterjtart of Eq. (27) after the finite element 
•li'cre! izai ion is the well-known stiffness expression 

A-.{r} = {P}-{f} (29) 

'.litii 'iie tansent stiffness tnatrix A’;, the vector of the incre- 
mental velocities {!'}. and the vector out-of-halance force rales, 
e.vter:iai force rates {P} minus internal force rates {f }. Theho- 
iiiocenou- ca-e of Eq. (29) indicates either the non-uniqueness 
of tin- e.juilihriuttt path at a stable or unstable bifurcation point 
or the iiiiujiie lull unstable situation at a limit point. Hence this 
criterion inav be evaluated by a determinant check or supple- 
mentarv eiuenvalue analysis for the load parameter parallel to 
tiie Ir.adill" process. 

E .< i, under the condition of static external loads and'slowly 
vrow iti” I r. ep . ffi.cis. the presence of snap-through buckling 
m.tki . tile inertia effects significant. In dynamic analyses, the 
Appli'-fl l.oilv forces include inertia forces. Assuming that the 
tna-- of the liodv considered is preserved, the mass matrix can 
he ev.'duated prior to the time integration using the initial con- 
niurat ii 'll. 


Finite rieinent solution of any boundary-value problem in- 
'.oi ve. !),e -ohition of the equilibrium equation (global) together 
with ilie constitutive equation (local). Both equations are solved 
diMuitane.'U.iv in a step by step manner. The incremental fonn 
of tlie stloi.al and local equations can be achieved by taking the 
intecratio!, over the incremental time step At = tj+i - tj. The 
fctantnilar rule has been apftlied to execute the resulting time 

iiiK-L’ration. 


C harlv. the numerical solution involves iteration. A simpli- 
ii*-d \>-rsion ol of Riks Wempner conslant-arc-length method 
has l,..pn utilized. This iteration procedure which is a general- 
izatiot; of the displacement control method also allows to trace 
the ri'.nlinear response beyond bifurcation points. In contrast 

0. the conventional Newlon-Raphson techniques, the iteration 
"f the iiietiiod takes place in the velocity and load rate spare, 
i io’ load step of the first solution in each increment is limited 
I ■■ coMfolling the length ds of the tangent. Either the length is 

1. e:u ron-tant in each step or it is adapted to the characteristics 
of 111. -olution. In each step the triattguiar-sized stiffness matrix 
i:a- to I..- checked for negative diagonal terms, indicating that a 
( rit ical p'liiii is reached. 


U . .ShaUow Circular Clamped Arch 

fh" theorv ami computational procedures described in the 
prec. dim sections have been applied to the creep buckling anal- 
of a -hallow circular clamped arch. The problem of the 
ciaiiiped arch beside being of some practical interest contains a 
:;u::ii'er of -iniilarities to that of the uniforndy loaded spherical 
Cap. file trend of results of the arch problem serves as a good 
liidicaior to the behavior of the latter. 


The shallow circular clamped arch subjected to a single cen- 
tral concentrated load, as shown in Fig. 1, is analyzed. The 
material chosen for the numerical experimentation is the car- 
bon steel C-45 (DIN 1720) with E = 10' psi, u = 0.3 and 
iTfi = 2.7 X 10^ psi at room temperature. The viscoplastic prop- 
erties (the scalar materia) functions) were obtained in Ref. 28. 



Fig. 1 CHamped Circular Arch 

The analysis is performed with the aid of 24 paralinear iso- 
parametric elements. Fig. 2. The paralinear isoparametric el- 
ement is intended for the analysis of oriented structures where 
the geometry is such that the thickness is small compared to 
other dimensions. The characteristics of the element are defined 
by the geometry and interpolation functions, which are linear 
in the thickness direction and parabolic in the longitudinal di- 
rection (see Fig. 2). (Consequently, the element allows for shear 
strain energy since normals to a mid surface before deformation 

remain straight, but not necessarily normal to the mid surface 
after deformation. 
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Fig. 2 Paralinear Isoparametric Element 

I 

The elastic behavior, corresponding to both axisymetric and 
t^ymetric response, is shown on Fig 3. These curves are in 
complete agreement with those produced byGjelsvik and Bod- 
ner [6], only because the young’s modulus and Poisson’s ratio 
values used are virtually the same (carbon steel C-45 here, and 
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2024-T4 aluminum alloy in Ref. 6 ). Note that these elastic 
response curves are hypothetical for our material but true for 
the 2024- T4 alloy. The true behavior for our material is elasto- 
' viscoplastic and it is labeled as such on Fig. 3. Note that this 
curve represents quasi-static (steady state) elastovisco- plastic 
response, as described by the chosen constitutive law. According 
to this, snapping occurs at a load of 26.20 lbs, primarily because 
of the low yield strength. Then, the post-limit point behavior 
seems to be primarily driven by viscoplastic behavior. 



Fig. 3 The Arch Response 

It is expected here that if loads up to approximately 14 lbs 
are reached quasi-statically and left applied for a long time the 
primary response will be creep and the critical time to creep 
will be extremely large. On the other hand, for loads between 
14 lbs and 26.2 lbs ( especially for the higher range ) the behav- 
ior will be a combination of creep and snap-through buckling. 
This is best demonstrated by the curve on Fig. 4. The applied 
load is reached quasi-statically in 13 minutes and then it is kept 
constant. The curve of Fig. 4 depicts the behavior of the arch 
in terms of midpoint deflection versus time. CVeep. initially, is 
very slow, then snap-through takes jilace in 32 minutes, curve 
IK', and then the creep behavior coiitiiiues until a critical time 
to creep (creep buckling occurs) is reached after a total time 
of 97 minutes. Note that for this loading condition the critical 
time to creep is 97 minutes. C'reep buckling and critical time to 
creep are based on the phenomenon that the deflection increased 
very rapidly. For loads higher than 26.2 lbs, it is expected that 
snapping will occur early, during quasi- static loading and then 
the creep behavior will be similar to that shown on Fig 4, past 
point c. 

The next example considers the influence of cyclic loading on 
the response. Figure 5 illustrates the load deflection behavior 
of the arch under a cyclically applied external load. The load 
is increased, quasi-statically, from zero to 26 lbs in 5 minutes, 
then it is held constant for 2.5 iiunutes, after that it is reduced 
to 20 lbs and held constant for 50 minutes, then raised Co 25.5 
lbs for 2.5 minutes and Anally it is reduced back to 20 lbs and 
held constant. 



Fig. 4 Deflection Time History 



Fig. 5 Multi-Cycle Arch Response 

The steady state response under this type of loading exhibits 
several limit points, which may imply that snapping is iimiii- 
nent shortly after the load reaches the value of 26 lbs (between 
points A and B on Fig 5). The dashed curve corresponds to 
the hypothetical elastic static response and it is only shown for 
comparison purposes. 

The last example presented in Fig 6. consider the influence 
of temperature on the arch behavior. The loading history is the 
same on the one shown on Fig. 4. The curve corresponding to 

T = 50"F' was discussed previously (Fig. 4), and it is used here 
as a basis for comparison. \Mten the temperature is raised to 
200"^ (after this, the loading is applied), thq time to snap is 
reduced to 26 minutes, while the critical time to creep is not 
appreciable affected. On the other hand at the highest tempera- 
ture T = 1000“F, for which results are shown. The critical time 
to creep is reduced to 62 minutes, and the steady state response 
does not show a clear snap-through behavior. 
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V. Discussion 

A' noted earlier, tlie main thrust of this work has been to - 
n* nioi.-trate i1,p effect of liighlv nonlinear material behavior on 
Uie -nap throiiah and creep buckling of shallow arches. It is ev- 
iii’ tit liia: In the presence of l)uth elastic and viscoplastic defor- 
ni.iiion til*' process of buckling assumes an entirely new charac- 
t-r. Hiirkling develops as a time- temperature dependent defor- 
m.iiion process under constant or variable loads of magnitudes 
»io;di< r ilian the elastic critical values. In arches under loads 
lo hm the critical values the structure initially deform quasi- 
'taiicaliv with the thentio- viscous terms manifesting themselves 
in the fortii of increasing displacetiieitt under, say a con.stant 
io,,.l. U'lien the ttiagttitude of the displacements reaches a cer- 
tain thre-hold state, the arch sna|>s dytiamirally into the post- 
lairklinu ronli"uration and then continues quasi-static defomia- 
tion a'.'aiti. 

The material constitutive relalioti has been proven to be ca- 
paide of reprodttcing the main characteristics of viscoplastic de- 
formatiotis. The ttiodiiied Riks/ AVettipner iteration scheme has 
I found to he a versatile technique in the pre-atid post criti- 

cal ranitc. 

Die influence of thermo- mechanical coupling can become 
verv large in stability problems. Such processes are always 
connected with a rapid growth of inhomogeneity of the state. 
Tliorough investigation of such problems, however, must he per- 
hirmeti with llie necessary detail. 
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